
###################
# IRT Fit Analyses
##################

# 1. Model fit -----------------------------------------------------------------
# Exploratory analysis:
modelfit1EA<-M2(mEA_yes_pos, type = "M2*", calcNULL = TRUE) # save statistics
modelfit2EA<-M2(mEA_yes_neg, type = "M2*", calcNULL = FALSE)
modelfit3EA<-M2(mEA_no_pos, type = "M2*", calcNULL = FALSE)
modelfit4EA<-M2(mEA_no_neg, type = "M2*", calcNULL = FALSE)
# Confirmatory analysis:
modelfit1CA<-M2(mCA_yes_pos, type = "M2*", calcNULL = FALSE)
modelfit2CA<-M2(mCA_yes_neg, type = "M2*", calcNULL = FALSE)
modelfit3CA<-M2(mCA_no_pos, type = "M2*", calcNULL = FALSE)
modelfit4CA<-M2(mCA_no_neg, type = "M2*", calcNULL = FALSE)

# Combine statistics into one data frame:
modelfit_list<- list(modelfit1EA, modelfit1CA,modelfit2EA,modelfit2CA, modelfit3EA, modelfit3CA, modelfit4EA,modelfit4CA)
model_fit <- bind_rows(modelfit_list, .id = "Model") # combine all statistics by model
model_fit<-model_fit%>% # rename models
  mutate(
  Model = case_when(Model == 1 ~ "IntPosEA",
                    Model == 2 ~ "IntNegEA", 
                    Model == 3 ~ "UnintPosEA", 
                    Model == 4 ~ "UnintNegEA", 
                    Model == 5 ~ "IntPosCA", 
                    Model == 6 ~ "IntNegCA",
                    Model == 7 ~ "UnintPosCA",
                    Model == 8 ~ "UnintNegCA",TRUE ~ NA))
row.names(model_fit) <- NULL # remove numbers next to name

# Create model fit table:
stargazer(model_fit, type= "latex", summary = F, 
           rownames = FALSE, digits = 3, digits.extra = 0) # table

# 2. Item fit ------------------------------------------------------------------
# Exploratory analysis:
itemfit1EA<-itemfit(mEA_yes_pos, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") # control for multiple testing
itemfit2EA<-itemfit(mEA_yes_neg, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 
itemfit3EA<-itemfit(mEA_no_pos, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 
itemfit4EA<-itemfit(mEA_no_neg, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 
# Confirmatory analysis:
itemfit1CA<-itemfit(mCA_yes_pos, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 
itemfit2CA<-itemfit(mCA_yes_neg, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 
itemfit3CA<-itemfit(mCA_no_pos, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 
itemfit4CA<-itemfit(mCA_no_neg, fit_stats = "S_X2", which.items = 1:8, p.adjust = "BH") 

# Combine into one data frame
itemfit_list<- list(itemfit1EA, itemfit2EA, itemfit3EA, itemfit4EA,
                    itemfit1CA, itemfit2CA, itemfit3CA, itemfit4CA)
item_fit <- bind_rows(itemfit_list, .id = "Model")

item_fit<-item_fit%>% # rename models
  mutate(
    Model = case_when(Model == 1 ~ "IntPosEA",
                     Model == 2 ~ "IntNegEA", 
                     Model == 3 ~ "UnintPosEA", 
                     Model == 4 ~ "UnintNegEA", 
                     Model == 5 ~ "IntPosCA", 
                     Model == 6 ~ "IntNegCA",
                     Model == 7 ~ "UnintPosCA",
                     Model == 8 ~ "UnintNegCA",TRUE ~ NA))
row.names(item_fit) <- NULL

# Table 1
item_fit_IntPos<-item_fit%>%
  subset(Model == "IntPosEA"|
           Model == "IntPosCA")
stargazer(item_fit_IntPos, type="latex", summary = F, rownames = FALSE, digits = 3, digits.extra = 0)

# Table 2
item_fit_IntNeg<-item_fit%>%
  subset(Model == "IntNegEA"|
           Model == "IntNegCA")
stargazer(item_fit_IntNeg, type="latex", summary = F, rownames = FALSE, digits = 3, digits.extra = 0)

# Table 3
item_fit_UnintPos<-item_fit%>%
  subset(Model == "UnintPosEA"|
           Model == "UnintPosCA")
stargazer(item_fit_UnintPos, type="latex", summary = F, rownames = FALSE, digits = 3, digits.extra = 0)

# Table 4
item_fit_UnintNeg<-item_fit%>%
  subset(Model == "UnintNegEA"|
           Model == "UnintNegCA")
stargazer(item_fit_UnintNeg, type="latex", summary = F, rownames = FALSE, digits = 3, digits.extra = 0)


# 3. IRT parameters ------------------------------------------------------------
# Exploratory analysis:
IntPosEA<-coef(mEA_yes_pos, IRTpars = TRUE, simplify = TRUE) # IRT parameters
IntNegEA<-coef(mEA_yes_neg, IRTpars = TRUE, simplify = TRUE) 
UnintPosEA<-coef(mEA_no_pos, IRTpars = TRUE, simplify = TRUE) 
UnintNegEA<-coef(mEA_no_neg, IRTpars = TRUE, simplify = TRUE) 
# Confirmatory analysis:
IntPosCA<-coef(mCA_yes_pos, IRTpars = TRUE, simplify = TRUE)
IntNegCA<-coef(mCA_yes_neg, IRTpars = TRUE, simplify = TRUE) 
UnintPosCA<-coef(mCA_no_pos, IRTpars = TRUE, simplify = TRUE) 
UnintNegCA<-coef(mCA_no_neg, IRTpars = TRUE, simplify = TRUE)

# Combine into one data frame:
IRT_list<- list(IntPosEA, IntNegEA, UnintPosEA, UnintNegEA,
                IntPosCA, IntNegCA, UnintPosCA, UnintNegCA)
IRT_fit <- bind_rows(IRT_list, .id = "items")

# Extract item statistics from each model in the list and bind them row-wise
IRT_fit <- map_df(IRT_list, ~{
  item_stats <- data.frame(t(unlist(.x$items)))
  colnames(item_stats) <- paste0("Item_", 1:ncol(item_stats))
  return(item_stats)
}, .id = "Model")

# Convert row names to a column
IRT_fit <- rownames_to_column(IRT_fit, var = "Parameter")

# Remove the numbers at the end of each row name
IRT_fit$Parameter <- gsub("\\.[0-9]+$", "", IRT_fit$Parameter)
IRT_fit$Parameter<- gsub("\\.", "", IRT_fit$Parameter)

IRT_fit <- IRT_fit %>% 
  dplyr::filter(Parameter == "a") 
  
IRT_fit <- IRT_fit %>% 
  mutate(
    Model = case_when(
      Model == 1 ~ "IntPosEA",
      Model == 2 ~ "IntNegEA",
      Model == 3 ~ "UnintPosEA",
      Model == 4 ~ "UnintNegEA",
      Model == 5 ~ "IntPosCA",
      Model == 6 ~ "IntNegCA",
      Model == 7 ~ "UnintPosCA",
      Model == 8 ~ "UnintNegCA", TRUE ~ NA_character_  # Ensure NA is character type
    ))

# Remove row names
row.names(IRT_fit) <- NULL

# Table:
stargazer(IRT_fit, type="latex", summary = F, rownames = FALSE, digits = 2, align = FALSE)

# 4. Person fit ----------------------------------------------------------------
# Exploratory analysis:
personfit1EA<-personfit(mEA_yes_pos)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2|Zh < -2)))
rownames(personfit1EA) <- c("Fit", "Misfit")
personfit1EA<-personfit1EA%>%
  mutate(Analysis="EA")
# Confirmatory analysis:
personfit1CA<-personfit(mCA_yes_pos)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit1CA) <- c("Fit", "Misfit")
personfit1CA<-personfit1CA%>%
  mutate(Analysis="CA")
# Exploratory analysis:
personfit2EA<-personfit(mEA_yes_neg)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit2EA) <- c("Fit", "Misfit")
personfit2EA<-personfit2EA%>%
  mutate(Analysis="EA")
# Confirmatory analysis:
personfit2CA<-personfit(mCA_yes_neg)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit2CA) <- c("Fit", "Misfit")
personfit2CA<-personfit2CA%>%
  mutate(Analysis="CA")
# Exploratory analysis:
personfit3EA<-personfit(mEA_no_pos)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit3EA) <- c("Fit", "Misfit")
personfit3EA<-personfit3EA%>%
  mutate(Analysis="EA")
# Confirmatory analysis:
personfit3CA<-personfit(mCA_no_pos)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit3CA) <- c("Fit", "Misfit")
personfit3CA<-personfit3CA%>%
  mutate(Analysis="CA")
# Exploratory analysis:
personfit4EA<-personfit(mEA_no_neg)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit4EA) <- c("Fit", "Misfit")
personfit4EA<-personfit4EA%>%
  mutate(Analysis="EA")
# Confirmatory analysis:
personfit4CA<-personfit(mCA_no_neg)%>% # check person fit using Zh 
  summarise(Zh.fit = prop.table(table(Zh > 2| Zh < -2)))
rownames(personfit4CA) <- c("Fit", "Misfit")
personfit4CA<-personfit4CA%>%
  mutate(Analysis="CA")

# Combine into one data frame:
personfit_list<- list(personfit1EA, personfit1CA, personfit2EA, personfit2CA,
                      personfit3EA, personfit3CA, personfit4EA, personfit4CA)
personfit <- bind_rows(personfit_list, .id = "Model")

# Convert row names to a column:
personfit <- rownames_to_column(personfit, var = "Fitting")

personfit<-personfit%>% # rename
  mutate(
    Model = case_when(Model == 1 ~ "IntPosEA",
                      Model == 2 ~ "IntNegEA", 
                      Model == 3 ~ "UnintPosEA", 
                      Model == 4 ~ "UnintNegEA", 
                      Model == 5 ~ "IntPosCA", 
                      Model == 6 ~ "IntNegCA",
                      Model == 7 ~ "UnintPosCA",
                      Model == 8 ~ "UnintNegCA",TRUE ~ NA))
row.names(personfit) <- NULL

# Remove the numbers at the end of each row name:
personfit$Fitting <- gsub("\\.[0-9]+$", "", personfit$Fitting )
personfit$Fitting <- gsub("\\.", "", personfit$Fitting)

# Table:
stargazer(personfit, type="latex", summary = F, rownames = FALSE, digits = 2)
